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Abstract 

An approach to derive low-complexity models describing thermal radiation for the sake 
of simulating the behavior of electric arcs in switchgear systems is presented. The idea is 
to approximate the (high dimensional) full-order equations, modeling the propagation of 
the radiated intensity in space, with a model of much lower dimension, whose parameters 
are identified by means of nonlinear system identification techniques. The low-order model 
preserves the main structural aspects of the full-order one, and its parameters can be straight¬ 
forwardly used in arc simulation tools based on computational fluid dynamics. In particular, 
the model parameters can be used together with the common approaches to resolve radia¬ 
tion in magnetohydrodynamic simulations, including the discrete-ordinate method, the P-N 
methods and photohydrodynamics. The proposed order reduction approach is able to sys¬ 
tematically compute the partitioning of the electromagnetic spectrum in frequency bands, 
and the related absorption coefficients, that yield the best matching with respect to the 
finely resolved absorption spectrum of the considered gaseous medium. It is shown how the 
problem’s structure can be exploited to improve the computational efficiency when solving 
the resulting nonlinear optimization problem. In addition to the order reduction approach 
and the related computational aspects, an analysis by means of Laplace transform is pre¬ 
sented, providing a justification to the use of very low orders in the reduction procedure as 
compared with the full-order model. Finally, comparisons between the full-order model and 
the reduced-order ones are presented. 


Keywords: Arc simulations. Radiative heat transfer, Model order reduction, Nonlinear estima¬ 
tion, Nonlinear model identification 

1 Introduction 

The switching performance of circuit breakers depends strongly on the behavior of the electric 
arc that originates when the contacts are opened in presence of relatively large electric current 
values [271E]. In turn, the arc dynamics are influenced by multiple interacting physical phe¬ 
nomena which, together with the short timescale of the arcing event and the large values of 
temperature and pressure, increase the complexity and difficulty of understanding, carrying out 
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experiments, and deriving numerical models of the switching behavior. Computational fluid dy¬ 
namic (CFD) approaches are being used in both public and private research efforts to simulate 
the time evolution of the plasma that carries the current during the interruption process, see e.g. 
miEiEiiia]. The CFD simulations are often coupled with solvers for the electro-magnetic (EM) 
phenomena, resulting in sophisticated multi-physics simulation tools (see mm) that allow one 
to gather an insight of what is actually happening during the current interruption process - 
aspects that are very difficult to quantify with direct measurements for the above-mentioned 
reasons. Such simulation tools provide a significant added value to explain the results of experi¬ 
mental tests and to support the development of switchgear devices, however they also bring forth 
an important issue in addition to the inherent difficulty of plasma modeling: the need to find 
a good balance between the accuracy of the employed physical models and their computational 
complexity. The modeling of radiative heat transfer during the arcing process is an illuminating 
example of such an issue. 

Radiation is one of the most important cooling mechanisms during switching, as it redis¬ 
tributes the heat produced by the current flowing through the plasma. Hence, accurate models 
of radiation are of fundamental importance to simulate the arc behavior, which is, due to the 
physics of radiation at the temperatures present in the plasma, a formidable task. Typically, 
the core of the arc is heated up to 25,000 K, meaning that the electromagnetic radiation emitted 
by ions, atoms, and molecules of several different species (e.g. nitrogen, oxygen, or copper) 
have to be taken into account. The relevant window of the electromagnetic spectrum ranges 
from 3 10^^ Hz-6 10^® Hz, corresponding to wavelengths between 10“® m and 5 10“® m, i.e. from 
infrared to ultraviolet. 

The main difficulty in simulating the electromagnetic radiation emitted by an arc derives from 
the complexity of the emission spectrum, where the relevant property, the absorption coeffi¬ 
cient, changes by many orders of magnitude at spectral lines of which several 10,000 exist in 
the range under consideration. The propagation of radiative heat in space for each frequency 
is modeled, under assumptions that are reasonable for the arcing phenomena encountered in 
switchgear devices, by a first-order differential equation taking into account the absorption and 
the emission of radiation along the direction of propagation. The energy removed from the arc is 
with this defined by the temperature, pressure, and composition distribution within. Due to the 
complexity of the emission spectrum, a simple discretization according to frequencies leads to 
hundred thousands of very thin frequency bands; within each one of such bands the absorption 
coefficient can be assumed to be constant for fixed temperature, pressure and composition of the 
gaseous medium. This however, leads to the same number of three dimensional field equations 
which needed to be solved. Given the large number of finite volumes that have to be considered 
in CFD simulations of realistic geometries (see e.g. [22])) the use of such a large-scale radiation 
model is not feasible. Hence, there is the need to derive models for the radiative heat transfer 
with much lower complexity, possibly without compromising too much the accuracy as com¬ 
pared with the large-scale model. This issue has been tackled by several contributions in the 
literature [Illl23l[I21EQ]. Most of the existing approaches consist in discretizing the fraction of 
the EM spectrum of interest into few bands, and assuming for each one of them some averaged 
absorption properties. The mentioned approaches have the advantage of being quite simple to 
implement, however in principle one should optimally choose ad-hoc different bands and aver¬ 
aged absorption coefficients as a function of pressure, temperature and chemical composition 
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of the considered medium, since the radiation parameters are affected by all these aspects. If 
the bands and/or the averaged absorption coefficients are not chosen in an appropriate way, 
the model accuracy is worse and more bands are needed to improve it, resulting in a relatively 
large number of bands (6-10) in order to achieve accurate results with respect to the original, 
large-scale model. Hence, this procedure can be time consuming, not trivial to carry out in a 
systematic way, and ultimately suboptimal in terms of complexity/accuracy compromise. 

In this paper, we present a new approach to derive small-scale, band-averaged models of 
the radiative heat transfer. We first describe the problem of radiation modeling from a novel 
perspective, where the aim is to approximate the input-output behavior of a large scale, linear- 
parameter-varying (LPV) dynamical system with that of a low-order one. The large scale 
system has one input (black-body intensity), one output (radiated intensity), three scheduling 
parameters (temperature, pressure and composition), and a large number of internal states (one 
for each considered frequency of the EM spectrum), while the low-order LPV system has the 
same input, output and scheduling parameters, but just a handful of internal states. From this 
point of view, the problem can be classified as a model-order reduction one m- Then, using 
classical tools for the analysis of signals and dynamical systems, we provide evidence that indeed 
models with quite low order (typically 2-3 bands) can be already good enough to capture the 
main behavior of the fnll-order model. Finally, we tackle the order reduction problem by using 
nonlinear system identification techniques (see e.g. m), where we collect input-output data 
from the large-scale system and use it to identify the parameters of the reduced-order model. The 
approach results in a nonlinear optimization problem (nonlinear program - NLP) with a smooth 
non-convex cost function and convex constraints, which are needed to preserve the physical 
consistency of the reduced-order model. We show through examples that the obtained reduced- 
order models enjoy a high accuracy with respect to the full-order one, while greatly reducing the 
computational times. As compared with the existing approaches, the method proposed here has 
the significant advantage of being systematic, i.e. there is no need to tailor or tune it for each 
different composition of the absorbing/emitting medium. The main user-defined parameter is 
the desired number of frequency bands in the reduced-order model, which can be then increased 
gradually until the desired tradeoff in terms of model quality vs. complexity is reached. 

The paper is organized as follows. Section [l . 1 1 provides a description of the problem we want 
to solve. In Section such a problem is analyzed from a system’s perspective and connec¬ 
tions are made to the order reduction of a large-scale LPV dynamical system. The proposed 
computational approach is described in Sectionfinally results are presented in Section]^ and 
conclusions and future developments are discussed in Section 

1.1 Problem formulation 

Let us consider a region in space containing a hot gaseous mixture of r G N different components. 
Each component is present in a fraction (e.g. of mass or of mole) G [0, 1], i = 1,... ,r,. We 
indicate with y = [yi,... ,yr]^ the composition vector of the medium. We consider a line of 
propagation along which we want to compute the radiated heat, and denote with x G M the 
position of a point lying on such a line. As the line crosses the hot matter, the temperature 
T(x) G M"*", pressure p{x) G M"*" and composition y(x) G [0, 1]'’ change with x (see Figure[^for a 
graphical visualization). The temperature, pressure and composition of the gas lie in some sets 
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Figure 1; Sketch of the considered problem. The aim is to compute the spatial distribution of the 
radiated intensity I{x,i') (dash-dotted line) for a given frequency v of the EM spectrum through a 
gaseous medium along a given direction with coordinate x. The temperature T (solid line), pressure p 
(dashed) and composition y (balloons of different colors) of the gas depend on x. The intensity at x = 0, 
is a known initial condition. 

of interest, T, V and y respectively, defined as follows: 

T = [Tmin, Tmax] c M+ 

'P — [Pmin) PmsLx] C M~*~ 

= \yG[o,iY-.f:yi 

i i=i 

Typical values defining the sets T and V are Tmin = 300 K, Tmax = 25, 000 K, p^m = 10^ Pa and 
Pmax = 10’^ Pa. 

Assuming without loss of generality that at x = 0 the radiated heat intensity, for a given 
frequency v of the EM spectrum, is equal to a given value Pfl, the distribution of the intensity 
as a function of x can be computed through the following ordinary differential equation: 

= a{T{x),p{x),y{x),iy){Ibb{Tix),u) - = Pfl, (2) 

where «(•,•,-, 1 ^) is the absorption coefficient and is the black body intensity, both 

pertaining to the frequency v. For a fixed value of a depends on temperature, pressure and 
composition, while Ibb is a function of temperature only. 

In equation ([^ , it is implicitly assumed that the radiation distribution has already converged 
to a steady state. Since the temporal dynamics of the radiated heat are much faster than the 
timescale of the phenomena of interest in arc simulations, such an assumption is reasonable. 
Similarly, scattering effects have been also neglected as they represent a negligible term for the 
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application considered here. For a more complete form of equation © the interested reader is 
referred to [26]. 

Thus, for each frequency v the corresponding ODE Q is characterized by two parameters, 
namely the black body intensity I^b and the absorption coefficient a. The dependency of the 
latter on T(x), p{x), y{s) renders the equation nonlinear. From a physical perspective, the black 
body intensity is a source of radiation, while the current intensity I{x,v) is a sink: the change 
of radiated intensity in an infinitesimal space interval dx is the difference between these two 
contributions, scaled by the absorption coefficient. The black-body intensity as a function of 
the frequency u and of temperature T is given by Planck’s law |26j : 

2 / 11/3 1 

Ibb{T,u) = —2- 

c _ I 

where h ~ 6.62 j g-i jg the Planck constant, c ~ 10® ms”^ ig the speed of light in vacuum 
and ks — 1.38 10“^® J ig the Boltzmann constant. The behavior of IbbiT,iy) as a function 
of V for some temperature values is shown in Figure]^ For the purpose of this study, the total 
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Figure 2: Planck’s law (j^ (solid lines) computed for different temperature values. The frequency band 
of visible light is also shown (between gray dashed lines). 

black-body intensity Ibb{T) per unit of surface and solid angle (black-body radiance) is also 
needed, given by Stefan-Boltzmann’s law: 

OO 

lbb{T) = I Ibb{T, v)du = (4) 

0 
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where usb — 5.6710“® Wm“^ K“^ is the Stefan-Boltzmann constant. From equation Q and 
Figureit can be noted that, for the temperature values of interest for arc simulations in circuit 
breakers (3,000K-25,000K), most of the black-body radiance (i.e. the area below the solid curves 
in Figure]^ is contributed by the frequencies in the range of approximately 3 10^^ Hz-6 10^® Hz. 

While the dependence of the black-body intensity on v and T is well-known from the above 
physical laws, the absorption coefficient a is a much more uncertain quantity. Gaseous media 
have absorption spectra (i.e. the function relating a to v) characterized by sharp lines at spe¬ 
cific frequencies, whose values depend on the composition of the mixture and whose number is 
typically very large. Moreover, the absorption coefficient at each of such frequencies depends 
strongly on temperature and (less markedly) on pressure. Broadening effects due to pressure 
are also important as they contribute to a spread of the absorption lines over the nearby fre¬ 
quencies. Physical models to capture such complex absorption spectra have been proposed in 
the literature, gum Eli. Moreover, databases of experimentally measured absorption data 
for several mixtures in specific frequency bands are available, typically for relatively low tem¬ 
peratures [2T]. In the following, we will assume that a is a known function of T, p, y, v, in the 
domain of interest; we call such information the “base data” T>-. 

V = {a{T,p,y,v)- VT eT,ype V, Vy G 3^, Vi/ G M+} (5) 

As an example of information contained in D, Figure shows the absorption spectrum for a 
mixture of 50% silver, 25% air and 25% hydrogen at 16,300 K and 10® Pa. For a fixed frequency 
of such a spectrum. Figure shows the corresponding absorption coefficient as a function of 
temperature and pressure (note the logarithmic scale for the absorption coefficient in this plot). 

Given the base data, the total radiated heat intensity Itot{x) is computed by integrating the 
intensity contributed by each frequency of the EM spectrum: 

OO 

Itot{x) = j I{x,v)di'. (6) 

0 

In order to obtain a tractable problem, a first step is to adopt a fine discretization of the 
frequency domain such that Q becomes a sum over a finite number of terms. Since the radiated 
intensity is relevant only in a specific region of the spectrum (compare Figure , the frequency 
discretization can be coarser outside such a region and finer inside, in particular around the peaks 


of the absorption spectrum. In this way, a finite number N of frequency values z/j, i = 1,..., 77 
are considered, each one being the middle point of an interval Az/j of the spectrum. Then, 
equation can be re-written as: 

N 

/(x) ~ ^/(x,z/i)Az/i. (7) 

i=0 

Equation Q, together with Q-Q and with the base data V ([^ evaluated at Vi,i = 1,..., A, 


form a high-dimensional model of the radiated heat intensity, named the Full-Order Model 
(FOM). The high dimensionality of the FOM comes from the fact that a large number N of 
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Figure 3: Example of absorption spectrum for a mixture of 50% silver, 25% air and 25% hydrogen at 
16,300 K and 10® Pa. 


frequencies is taken into account in the discretization, so that the approximation error is small. 
Typical values of N for the conditions of our interest are in the order of 1-2 10®. Due to the large 
number of considered frequencies, the FOM can be used effectively only in very simple cases, 
for example in one-dimensional problems. In fact, in order to solve the radiative distribution in 
two- and three-dimensional cases, as it is needed for example in CFD simulations of real devices, 
one would have to solve many equations (whose nature depends on the chosen method, like the 
so-called P-N methods, the discrete ordinate method or the photohydrodynamic approach, see 
|26j . [To]) for each one of the considered frequencies, leading quickly to an intractable problem. 
In this paper, we tackle this issue by deriving a method to compute low complexity models for 
the radiative heat transfer. More specifically, we consider the following problem: 

Problem 1 ; Given the full-order model defined by ©-(HI), 0 and the base data V, derive 
a model with the same structure, i.e. where the total intensity is the sum of a finite number 
of contributions obtained by partitioning the frequency domain, but where the number of such 
partitions is very small, while still capturing accurately the total radiated heat intensity. 

We call such a simplified model the Reduced Order Model (ROM). The ROM can then be 
effectively used to model the propagation of heat via radiation in many applications, including 
full three-dimensional simulations of plasma arcs encountered in switchgear devices. 

We remark that, in light of Problem]^ in this work we will consider the radiated intensity 
computed with the FOM as “exact”. We will thus evaluate the quality of a given ROM by 
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Figure 4: Example of absorption coefficient for a mixture of 50% silver, 25% air and 25% hydrogen as 
a function of temperature and pressure, for v = 4.57110^^ Hz. 


assessing the discrepancy between the radiative heat intensity given by the latter and the one 
given by the FOM. In other words, we don’t consider here the accuracy of the FOM with respect 
to the real-world behavior of gaseous media. Indeed, the topic of modeling accurately and/or 
measuring the absorption spectrum of a given gas as a function of temperature and pressure 
(i.e. to compute the base data 2?) is by itself an important and active research area [41IT61 l3l IH] . 
however it is outside the scope of this work, which is focused on the simplification of the FOM 
into a computationally tractable model. On the other hand, the method presented in this paper 
does not depend on the specific absorption spectrum, i.e. it can be applied systematically to any 
base data, and it yields quite small discrepancies between the derived ROM and the employed 
FOM, such that the error between the ROM and the real behavior of the considered medium 
depends ultimately only on the quality of the FOM. 

2 A system’s perspective of radiative heat transfer 

2.1 Equivalent input-output models of the radiative heat transfer 

As a preliminary step to address Problem we re-write the FOM in a slightly different form, 
which is convenient to show that this model can be seen as a single-input, single-output (SISO) 
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Linear Parameter-Varying dynamical system. Let us define the spectral emissivity e(T, u) as: 


e(r, u) 


hb{T, v) 
~hb{T) 


( 8 ) 


i.e. the ratio between the total black body intensity and the one pertaining to each frequency 
of the EM spectrum. By construction we have e(T, v) G (0,1), V(r, u) and e(T, u)di> = 1 VP. 
The course of e(T, v) is illustrated in Figure 

By inserting the spectral emissivity Q in equation Q, the radiated heat intensity for each 
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Figure 5: Spectral emissivity e(T, v) as a function of EM frequency (solid lines), computed for different 
temperature values. The frequency band of visible light is also shown (between gray dashed lines). 

frequency value Vi,i = 1,... ,N of the FOM can be equivalently computed as: 

= (^{T{x),p{x),y{x),Vi) {e{T{x),i'i)lbbiT{x)) - I{x,i'i)) , 1(0, z^*) = (9) 

Equation Q is just a re-writing of Q but, together with equation Q, it highlights the fact 
that the FOM can be seen as a dynamical system where x is the independent variable, T{x) is 
an exogenous scalar input, /(x, Vi),i = 1,..., V are N internal states, Itot{x) a scalar output. 
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and p{x),y{x) are space-dependent parameters. In virtue of the fact that the function relating 
the temperature T to the total black-body intensity Ihb is known, one can consider the latter as 
input to the system hence putting into evidence the Linear Parameter Varying (LPV) structnre 
of the model: 


= A{T{x),p{x)Mx))I{x) + B{T{x),p{x)M^))hmx)) 

hot = Cl{x), 

where 

l{x) = [I{x,ui),...,I{x,UN)f 
MT,P,y) = diag(^[-a{T,p,y,ui),... ,-a{T,p,y,UN)f^ 

B{T,p,y) = [a(r,p,y,i/i)e(r, j/i),...,a(r,p,y,z/Ar)e(T,i/Ar)]'^ 


( 10 ) 


( 11 ) 


A block-diagram of the FOM from this new point of view is shown in 


Figure [^a). 


Such a 


4f)(i 


dl{x, Vil 
dx 


a(7’,p,y,Vi)(e(7’,Vi)4f, -/(x,Vil ) 


Avi 


• N 


IfOt 


dI(x,VN] 

dx 


aiT,p,y,VN)(_e(T,VN)hb- I(x,Vn] ) 


!\v 


N 


4f)(i 


dl(x, pi! 
dx 


= d(T,p,y,p.i\ (e(F,p,y,Pi)4h -/(x,pil ) 


• M «N 




^tot 


di{x, Pm\ 
dx 


a\ 


T, p, y, PmI (e(Fi Pi y, PM) 4 t) - Pm\ ) 


Figure 6: Equivalent block diagram of the (a) full-order model, FOM, and (b) reduced-order model, 
ROM. The dependence of T,p,y on x is omitted for the sake of readability. 
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system’s perspective of the radiative heat transfer equations is the fundamental step at the basis 
of all the developments described in the remainder of this paper. We can now introduce more 
precisely the structure of a candidate reduced-order model (ROM) meant to approximate the 
FOM. Let us denote the bands of the ROM with ^i, ... ,hmi with M N (e.g. M = 2), and 
the corresponding vector of inner states with I{x) = [I(x,/ii),... Then, 

we can write the equations describing the ROM as: 


dl{x) 

dx 


= MTix),pix),y{x)) I{x) + B{T{x),p{x),y{x)) hb{T{x)) 


( 12 ) 


Itot{x) = c I{x), 


where 


MT,P,y) = diag {[-a{T,p,y, pi),...,-a{T,p,y, pn)] J 

B{T, p, y) = [d(T, p, y, pi)e(T, p,y,pi),..., d{T, p, y, pn)KT, P, y, G 

C = [1,.. ., 1] G 


(13) 


plxM 


Figurej^b) gives a graphical representation of the ROM. Practically speaking, from (12)-(13) one 
can see that each component of vector I(x) accounts for a certain portion of the total radiated 


heat intensity, in complete analogy with the FOM (10), the only difference being the number 


of frequency bands, which in the ROM is much smaller than in the FOM. For each band pi, 
the parameters d{T,p,y, pi) and e{T,p,y, pi) have thus the meaning of “equivalent” absorption 
coefficient and spectral emissivity, respectively. The task of deriving a ROM for the radiated heat 
transfer from the FOM is referred to as model order reduction. Computing a ROM is equivalent 
to assigning suitable values to d and e as a function of the underlying parameters T,p and y. In 
particular the collection of all bands has to form a non-overlapping partition covering the whole 
EM spectrum. Then, the value of e{T,p,y, pi) has to correspond to the integral of the spectral 
black body intensity ([^ over the frequency band pertaining to pi. Hence, one can equivalently 
state that computing a ROM amounts to choose, for each pair of pressure and composition 
values, a partition of the EM spectrum (which defines the equivalent emissivity as a function 
of temperature) and the courses of the corresponding equivalent absorption coefficients as a 
function of temperature. In previous contributions in the literature, e.g. |2n[ I23| . the task of 
defining the ROM has been carried out by picking a finite number of bands covering the EM 
spectrum and then computing the equivalent absorption coefficients d{T,p,y, pi) through some 
averaging procedure on the portion of the absorption spectrum contained in each band. This 
approach is simple to implement but it has the drawback of not being systematic, since both the 
choice of the band cuts and the averaging of the absorption coefficients have to be made by the 
user, without an immediate link to the accuracy of the resulting ROM. In the next sections, we 
will present a new order reduction approach to compute the partitioning of the EM spectrum 
and the values of d{T,p, y, p) in a systematic way, that yields quite accurate results as compared 
with the FOM. 

About this last point, i.e. the accuracy of the ROM, we note that the input of the FOM 
and of the ROM is exactly the same, corresponding to the total black-body radiation for the 
considered temperature profile, Ibb{T{x)). Therefore, it is quite intuitive that the discrepancy 
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between the total intensity given by the FOM, Itot{x)-, and the one predicted by the ROM, 
Itot{x), for the same temperature profile T{x) (i.e. the same distribution of Ibb{T{x))), represents 
a reasonable indicator of the accuracy of the reduced order model. In other words, the error 
signal A/(x) = Itot{x) — Itot{x) will be considered to evaluate the goodness of a given ROM. 
This choice is motivated by the fact that the total radiated heat intensity is the main quantity of 
interest predicted by the ROM when it is embedded in multi-physics simulations of arc plasma, 
since it is used to compute the heat transferred from the plasma volume to the walls, and also 
(through its divergence) the heat redistributed within the plasma volume. Hence, the ROM 
should reproduce this quantity as accurately as possible with respect to the FOM, given the 
same spatial distribution of temperature, pressure and chemical composition. 

Before going to the details of the proposed method to derive the ROM, a sensible question 
to be addressed is whether the approximation problem we are dealing with has a reasonably 
good solution or not. More specifically, recall that we aim to approximate the input-output 
behavior of a large scale system, with hundred of thousands of internal states, with that of a 
small-scale one, with at most a handful of internal states. It is not immediately clear if there 
exist such a low-order ROM still capable of delivering high approximation accuracy, since this 
aspect depends on the characteristics of the FOM, i.e. on the underlying physics of the radiative 
heat transfer. In the next section, we exploit the system’s perspective described above to provide 
an intuition that indeed a ROM with a handful of bands can capture most of the input-output 
behavior of the FOM. 


2.2 Frequency-domain analysis 


Analyzing the complexity of the FOM in its form (10) in the domain of position x is not 


straightforward due to the very large number of internal states, each one following its own 
dynamic evolution. Besides noticing that the FOM is given by the sum of a large number 
of non-interacting, asymptotically stable first-order systems, all driven by the same input and 
whose (position dependent) poles are given by —a{T,p,y), little less can be said. 

However, if we consider fixed values of temperature, pressure and composition, T,p and 
y respectively, and we assume that only infinitesimal perturbations of temperature take place 
in space, such that the absorption properties of the medium can be assumed constant, we 
immediately notice that the FOM becomes a linear-parameter-invariant (LPI) system, to which 
well-assessed tools in systems theory and signal processing can be applied. In particular, after 
establishing the analogy between the position x in the FOM with the continuous time variable 
in dynamical systems, we can study the input-output response of the FOM to such inhnitesimal 
temperature variations by applying the Laplace transform [25] to its equations and deriving the 
transfer function G{s) from its input Ibb{s) to its output Itot{s), where s is the Laplace variable: 


G(s) = |^ = C(sI-A(r,p,y)) 'R(r,p,y), 


(14) 


In (14), A^B are the matrices given in equation 0 and evaluated at the chosen temperature, 
pressure and composition values, I is the identity matrix of suitable order and (si — A)~^ denotes 
a matrix inverse operation. A tool commonly used to analyze the transfer function of a dynamical 
system is the Bode diagram of the corresponding frequency response, obtained by evaluating the 
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magnitude and phase of G{juj), where ui = ^ assumes here the physical meaning of the frequency 
of purely sinusoidal oscillations (in space) of the black-body intensity (i.e. of temperature), with 
infinitesimal amplitude, with the period r measured in m (e.g. u) = 628 rad/m corresponds to a 
period of oscillation of the input of roughly 10“^ m). As an example, the Bode diagram of the 
FOM frequency response obtained by fixing T = 16, 300 K, p = 5 10^ Pa and y containing 50% 
silver, 25% air and 25% hydrogen (whose absorption spectrum is shown in Figure is shown 
in Figure It can be noted that the overall behavior of the FOM is that of a low-pass filter, 




Figure 7: Example of frequency response for a mixture of 50% silver, 25% air and 25% hydrogen at 
16,300 K and 10® Pa. Solid: full-order model; dashed: reduced-order model with M = 5 bands. 

whose frequency response is shaped by the contributions of the large number of internal states 
present in the system. As a quantitative example regarding this case, a sinusoidal oscillation 
of the black-body intensity (due to an oscillation in temperature) with a period of 27r 10^ m 
would give rise, in the gaseous medium, to a sinusoidal distribution of radiated heat intensity 
with unchanged amplitude and phase with respect to what would happen in vacuum (compare 
Figure [^with w = 10“^ rad/m), i.e. the amplitude of the oscillations would be equal to that of 


13 







the input black-body intensity, and the spatial distribution would show almost zero phase shift. 
On the other hand, the spatial distribution of radiated heat induced by a sinusoidal oscillation 
of the black-body intensity with a period of 27r m (compare again Figure]^ with a; = 1 rad/m) 
would have an amplitude equal to only about 23% with respect to that of the input, with a 
phase lag of about 45 deg. 

If we consider now the order reduction of such a dynamical system, we see that this is a 
standard problem of model order reduction of LPI systems, for which a well-established literature 
exist, see e.g. Thus, we can use one of the existing approaches to derive a reduced order 

model that approximates the FOM for the chosen values of T,p and y. After deriving the ROM, 
the related transfer function can be computed as (compare equation (HI): 

Gis) = i^ = c(sI-AiT,p,y)Y'B{T,p,y). (15) 

hb{s) ^ ^ 

A comparison between the Bode diagram of the latter and that of the FOM reveals that up 
to very small gain values of about such that the corresponding radiated heat intensity 

is negligible, the frequency responses of the two models are practically super-imposed, hence 
showing a very good agreement between them. Moreover, applying this procedure for many 
values of T,p and y, chosen by gridding their respective domains T, V and T, shows that such 
a good agreement is obtained always with no more than four-five bands in the ROM. A good 
agreement up to a gain of about 5% is obtained with just two bands in the ROM, for oscillation 
periods of fractions of millimeters. Indeed, this level of accuracy would be enough for the 
sake of arc plasma CFD simulations, where the resolution of the spatial discretization of the 
considered volume is of the order of 10“^ m. An example of the obtained results is depicted in 
Figure too. In particular, the ROM whose frequency response is shown in the hgure has five 
bands, pi,..., and the corresponding absorption coefficients and emissivities (and frequency 
boundaries) are equal to 


a(r,p,y,Mi) 

= l.llO-^m-i 

e(T,P,y,lii) 

= 7.2 IQ-i 

[0, 1.6 1033] Hz 


a{T,p,y, P 2 ) 

= 9.210-im-i 

KT,P,y, P 2 ) 

= 2.1 IQ-i 

(1.61033, 2.41033] 

Hz 

&{T,p,y, ps) 

= 1.310^m-i 

HT,p,y, ps) 

= 5.4 10-2 

(2.41033, 3.1 1033] 

Hz 

a{T,p,y, Pa) 

= 1.2102m-i 

e{T,p,y, Pa) 

= 1.2 10-2 

(3.1 1033, 3.8 1033] 

Hz 

a{T,p,y, pA 

= l.TlO^m-i 

KT,P,y, Pb) 

= 3.9 10-3 

(3.8 1033, -hoo)Hz 



A comparison between these values and Figure shows that the values of the ROM absorption 
coefficients correspond to the dominant poles of the FOM, as the intuition would suggest. 

Overall, the analysis reported so far provides if not a rigorous proof at least an indication 
that the problem we are dealing with has a reasonable solution, using ROMs of quite low order. 
In section we describe in details the solution approach that we propose to deal with the 
linear-parameter-varying case. 

2.3 Position discretization 

Before proceeding further, it is convenient to introduce the discretized versions of the FOM and 
of the ROM, where the position variable x is taken at nodes Z G N, that are equally spaced 
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by an interval Ax. The latter has to be chosen according to the features of the problem at 
hand, trading off computational speed with a sufficiently fine discretization, which can capture 
well the fastest transients of the model’s input and output. As a rule of thumb, Ax ~ (5/10 can 
be chosen, where <5 is the smallest space-scale of interest in the problem (typically in our case 
6 ~ 10“^ m, as mentioned above). Hence, all the space-dependent variables (i.e. T,p, y) are now 
evaluated at discrete position values. The discretization of the radiation models is carried out 
by assuming that such variables are constant between two subsequent position nodes, xi and 


xi+i = xi + Ax, and then computing the explicit integration of (10) (for the FOM) and (12) 
(for the ROM). In particular, for the full-order model we have: 


I{xi+i) = Adp"{xi),p{xi),y{xi)) I{xi) + Bd{T{xi),p{xi),y{xi)) Ibb{T{xi)) 
Itot{xi) = C I{xi) 


where the matrices Ad and Bd are computed as: 

^d{T,p,y) = diag([a{T,p,y,iyi),...,a{T,p,y,UN)Y 


t,NxN 


Bd{T,p,y) = [{1 - a{T,p,y,i^i))e{T,i^i),... ,{1 - a{T,p,y,i^N))eiT,i^N)] G 


pNxl 


(16) 


(17) 


and 


a{T,p,y,Vi) = e 


— p-a(T,p,y,Vi)/\x ■ _ 


, Z = 1, . . . ,M. 


The matrix C in (16) is the same as in (10), since the output equation is static and thus it is 


not changed by the discretization of the position x. Similarly, for the reduced-order model we 
have: 


I(x/+i) = Ad{T{xi),p{xi),y{xi)) I{xi) + Bd{T{xi),p{xi),y{xi)) Ibb{T{xi)) 

itot{xi) = Cl{xi) 


(18) 


where the matrices Ad and Bd are computed as: 

Ad{T,p,y) = diag([a{T,p,y,pi),...,a{T,p,y,pM)f 


G 


pMxM 


Bd{T,p,y) = [(1 - a{T,p,y, pi))e{T,p,y, pi), ..., (1 - a{T,p,y, pM))e{T,p,y, pM)f G 

(19) 

( 20 ) 


and 




Also for the ROM the matrix C in (18) is the same as in (12). 


The reason why we employ such a space discretization is twofold: on the one hand, it is 
needed to obtain a finite-dimensional computational problem, on the other hand it improves 
the computational efficiency of the method (in particular by using a constant discretization step 
Ax). 


3 Model order reduction of the radiative heat transfer equation 

3.1 Solution approach 


Considering the analysis of Section 2.2, one can be tempted to derive the ROM, i.e. the functions 
d{T,p,y) and e{T,p,y), by gridding the domains T,V and V and for each triplet {T,p,y) 
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compute the equivalent absorption coefficient and emissivity of the corresponding LPI model, 
using well-assessed and efficient model order reduction techniques. Then, the ROM can be 
obtained by interpolating among the computed reduced-order LPI models. This approach could 
work well if only pressure and composition dependence were considered, since the sensitivity of 
the base data on these values is mild and hence one can be confident that interpolating among 
the bands computed at different triplets yields correct results. In other words, for a given 
temperature T, the absorption spectrum of the FOM does not change dramatically between 
two neighboring pairs (pi,yi) and {P 2 ,y 2 ) within the pressure and composition intervals of 
interest for the application considered here, so that for each band pi of the ROM it is safe to 
interpolate between the coefficients a{T,pi,yi, fii) and a{T,p 2 ,y 2 , pi)-, computed independently 
by means of LPI model order reduction, to obtain the ROM coefficients for generic values 
of {T,p,y) with p G [Pi; P 2 ]) V ^ [Viy ^ 2 ]' However, this procedure would not achieve good 
results when temperature dependence is taken into account as well. In fact, the dependence of 
the absorption spectrum on temperature is very strong (compare Figure Q, so that the ROM 
coefficients pertaining to the same band (e.g. pi) but computed at two different temperature 
values, even with the same pressure and composition values, might be completely unrelated 
to each other, and interpolating between them can give highly inaccurate results. Driven by 
these considerations, we adopt a hybrid strategy, where we grid the domains V and y and 
for each pair (p, y) we derive the partition of the EM spectrum in M bands and the related 
functions d{T,p,y,pi), i = 1,...,M that define the ROM. Since the corresponding FOM is 
now parameter varying (because we let the temperature distribution change while fixing only 
pressure and composition), this problem falls in the class of model-order reduction of LPV 
systems, for which, differently from the LPI case, few results exist in the literature [281124j . and 
their practical applicability to systems with ~ 10^ states, like the FOM in our problem, is not 
straightforward. For the above reasoning, in the remainder of this section it is assumed, unless 
otherwise stated, that a fixed pair (p, y) of pressure and composition values has been chosen, 
and that the only space-varying variable is the temperature T. The complete ROM can be then 
obtained by repeating the LPV order-reduction for all the pairs (p, y) chosen by gridding the 
respective domains, and then interpolating among the obtained values to compute the equivalent 
absorption coefficients and frequency bands for a generic triplet (T,p, y). 

Remark 1 The simplification introduced by fixing pressure and composition is made possible 
thanks to the above-discussed particular properties of the absorption spectra of the considered 
gaseous media. Indeed, such a simplification by itself is not required for the order-reduction 
technique described in the following, which may straightforwardly be applied also with varying 
pressure and composition, hut at the price of higher computational requirements. In our experi¬ 
ence, using constant pressure and composition in the order reduction computation yields accurate 
enough results for the applications of interest. 

In the following, we propose to address the LPV order reduction problem with a nonlinear 
system identihcation approach (see e.g. m), in which we search, within a given set of possible 
ROMs, the one which is closest to the FOM according to a pre-defined optimality criterion. 
In the next sections, this task is brought in the form of a tractable optimization program. It 
has to be noted at this point that, in the literature on system identification, there exist several 
contributions devoted to the problem of identification of LPV systems, see e.g. [5] and the 
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references therein. However, such results are not applicable in our case, due to the additional 
constraints that are present on the ROM, namely the need to preserve a specific structure 
where the total radiated intensity is the sum of the contributions given by the frequency bands. 
Differently from such previous approaches, the one proposed here is able to take into account 
these constraints, since it is specifically tailored for the considered application. 


3.1.1 Cost function and model set 


In order to have an optimization problem that can be solved with common numerical techniques. 


two main ingredients need to be defined: the set of reduced-order models of the form (18) where 


we carry out our search, denoted with % (model set), and a cost function J giving a measure 
of how much a given ROM H £ T-L is close to the FOM. The model set Ti should represent 
the limitations that we want to impose on the ROM in order to account for the physics of 
the problem. A ROM is fully characterized by its parameters a{T,p,y, Hi), e{T,p,y, pi), i = 
1,..., M, which in the considered settings are, for each band pi, functions of temperature only. 
To be consistent with the underlying physical phenomena, the equivalent absorption coefficients 
should be positive (to retain an asymptotically stable model): 


a{T,p,y,pi) > 0, VT E T, i = 1,... ,M 


( 21 ) 


and the equivalent emissivities should lie in the interval [0, 1] and sum to one over all the 
considered temperature range, in analogy with (Isl), so that the black-body limit is not violated: 


0 < e{T,p,y,pi) < 1, VT E T, i = 1,. 




M 


E KT,P,V,Pi) = 1, VT E T 


( 22 ) 


i=l 


Since equation (20) is invertible, we can equivalently consider the functions a{T,p, y, p) to define 
the model set. We select this alternative for the sake of computational efficiency, as we will 


discuss more in details in section 3.2.3 Then, the constraints (21) can be re-written as: 


0 < a{T,p,y,pi) < 1, VT E T, i = 1,... ,M 


(23) 


As a final step to define T-L, we choose a finite parametrization of functions a{T,p,y, pi) and 
e{T,p,y, Pi), in order to obtain a hnite dimensional model set (hence also a finite dimensional 
optimization problem): 


a{T,p,y,p) = faiT,0a{p,y,p)), 0a{p,y,p) E 
e{T,p,y,p) = fe{T,ee{p,y,p)), ee{p,y,p) E 


(24) 


where fa, fe are chosen, once again, to tradeoff the flexibility of the parametrization with com¬ 
putational efficiency. In particular, a convenient choice for /„ is the class of piecewise affine 
functions of temperature, while fg is taken such that the corresponding parameters Og define 
the partitioning of the EM spectrum in a finite number M of bands. These choices have the 


advantage of yielding a convex model set (see section 3.2.2 for details), hence improving the 
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computational efficiency and stability of the procedure. For a given pair {p, y), let us collect the 
parameters 9a and Og for all the ROM bands pi in a single vector 9: 


9 = 


^dip,y, pm) 

&e{p,y,pi) 


(25) 


OeiP,y,PM) 

where ne = M{n 0 ^ + n$.) is the total number of optimization variables. Then, considering 


(22)-(25), we can define the model set % as 


n = 


9 £ 


0 < faiT, 9g{p, y, Pi)) = 

0 < fe{T, 9g{p, y, Pi)) < 1, VT e T, i = 1,..., M 


M 


VT e r 


E fe{T,9e{p,y,Pi)) = 1 

i=l 


(26) 


As to the cost function J, this should represent a criterion by which the accuracy of the ROM 
is evaluated. As discussed in section an indicator of the accuracy of a ROM, suitable for the 
considered application, is related to the error between its total output intensity profile and that 
of the FOM, given the same temperature profile. Motivated by this consideration, we select 
a series of temperature profiles, and then evaluate the discrepancy between the corresponding 
outputs generated by the FOM and those given by the ROM. More specifically, let us consider 
a space interval X = [0,x], chosen such that x = (L — 1) Ax for some L G N. Similarly to 
the choice of Ax, the value of x depends on the problem at hand: a typical choice is three-four 


times the largest space-scale of interest (see section 3.2.1 for details). Thus, in the interval X we 
have L position nodes xi, I = 1,..., L, with xi = 0 and xl = x. Consider now a finite sequence 


of Lt G N discretized temperature profiles Tj = [Tj{xi),... ,Tj{xL)\ G 




, j = 1,.. .,Lt. 


In face of each of such temperature profiles, the FOM provides a corresponding profile of the 
radiated intensity, itot,j{Tj) = [Itot,j{xi),... ,Itot,j{xL)\^ G j = 1,...,Lt, and similarly 

the ROM, for a given value of 9 provides an approximated one, Itot,j{Tj,9). Let us dehne the 
weighted error profile as 


AIj{9) — Wj{Itot,j{fj) itot,jiTj, 9)), j — 1,, Lt 


(27) 


where wj > 0 is a scalar weight (the specific choice of Wj is discussed in section 3.2.1). Then, 
we take the cost function as: 

L'j' 

J{9) = ^AlJ{9)AI,{9), (28) 

i=i 


i.e. the sum, over all error profiles and over all position nodes of each profile, of the weighted 
intensity error squared. 
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3.1.2 Nonlinear Program formulation 

We can now write the model order reduction problem in a computationally tractable form: 


min J(0) 
e&H 


(29) 


where T-L and J{6) are given in (26) and (28), respectively. Namely, the aim of the optimization 


problem (29) is to search, within the model set the value of the parameter vector 6 that 


achieves the best fitting criterion J{6). Problem (29) is a finite-dimensional nonlinear program 


(NLP) which, for a suitable choice of the model parametrization (24), can be tackled with state- 
of-the-art numerical methods m- As we show in the examples of section a typical dimension 
of the optimization variable is ng ~ 77, corresponding for example to a 3 -order model with 25 
parameters defining the functions a{T,p,y, fii) and two parameters defining the frequency cuts 


in the EM spectrum that separate the three bands pi, p, 2 , ^ 3 - In general, the cost function (28) 


is non-convex with respect to the decision variable 9, so that only a local solution to problem 


(29) can be computed efficiently with such problem dimensions. As we show in the examples of 
section 0. the obtained solutions nevertheless yield ROMs that are very accurate with respect 
to the FOM. 


In the next section, we discuss several aspects involved in the formulation and solution of (29) 


including the design of the input temperature profiles and the choice of the weights in the cost 


function J (28), the choice of model parametrization, the numerical approach to solve the NLP, 


finally some ways to improve the efficiency and stability of the numerical optimization. 


3.2 Computational aspects 

3.2.1 Input design, weighting coefficients and initial conditions 

From an application’s perspective, the position interval x, the resolution of the position dis¬ 
cretization Ax, and the temperature profiles Tj, j = 1,..., Lt should be representative of the 
typical temperature distributions in the plasma that is generated during the switching process. 


About the choice of Ax, considering the discussion in section 2.3, a reasonable tradeoff between 
accuracy and computational speed is Ax = m, while for x a good choice in our experience is 
X = 2 10“^ m, for arcs whose width is between ~ 3 10“^ m and ~ 2 10“^ m, corresponding to cur¬ 
rent values in the range 1-40 10^ A. Finally, for the temperature profiles, we select a bell-shaped 
function, where we can adjust the maximum temperature reached as well as the steepness of the 
rising and falling slopes. Then, we generate a series of such profiles by cycling through differ¬ 
ent maximum temperatures and transient slopes. Furthermore, from a system identification’s 
perspective, given the nonlinearity of the system’s equations, different initial conditions for the 
internal states of the model should be considered, and the input should be designed in order 
to excite the system in a broad range of frequencies. For the former aspect, we replicate the 
temperature profile several times, so that the system is presented more than once with the same 
temperature profile, but each time starting from the internal states resulting from the previous 
profile. For the second aspect, white noise processes are a well-known choice m to excite the 
system’s dynamics over all frequencies, hence we super-impose to the computed profiles a uni¬ 
formly distributed white noise temperature signal whose amplitude is a fraction (e.g. 25%) of 
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the highest temperature in the original profile. With these choices, typical temperature prohles 
are shown in Figure [^a)-(b), together with the resulting intensity distributions computed with 
the FOM for a mixture of pure air. 


(a) 


(b) 




Position X (m) Position x (m) 


Figure 8: Example of input temperature profiles used in the order reduction procedure (solid gray lines) 
and radiated heat intensity distribution computed with the full-order model (dashed black lines). The 
maximum temperature values prior to adding white noise are (a) 10^ K and (b) 1.7 10^ K . Base data: 
pure air at 1.5 10® Pa. 


Due to the fourth-order dependence of the total black-body intensity on temperature (see 
il»: the radiated intensity obtained with different temperature profiles can have different orders 
of magnitude, compare for example Figure [^a) with Figure [^b), where the maximum temper¬ 
atures are about 10^ K and 1.7 10^ K, respectively. If no corrective measure is taken, such a 
difference among the various intensity profiles would lead to poor accuracy of the reduced order 
models, due to the biasing towards the higher intensity values in the cost function. In order 
to compensate this effect, we select the weights Wj in (27) on the basis of the radiated heat 
intensity given by the FOM when the corresponding temperature prohle, T,-, is considered. In 


this way, all the error profiles (27) are normalized with respect to the corresponding intensity; 
typical specific choices for the weights include the maximum or the average radiated heat, i.e.: 

Wj = max(Itotj), j = I,..., Lt 
or 

Wj = mean{Itot,j), j = 1 ,...,Lt 

Finally, it is worth mentioning the choice of the initial conditions for the radiated heat 
intensity, for both the FOM and the ROM, i.e. the values of vectors I{xi) and I{xi) needed 


to simulate the models (16) and (18), respectively, in the interval X = [0, x]. To this end, we 
assume that for x < 0 (i.e. outside such an interval) the incoming radiation in positive x direction 
corresponds to the steady-state intensity at ambient temperature Tq, i.e. Ibb{Ta, i = , X 


20 


Radiated heat intensity (W m"'^ rad’ ) 













































for the FOM and 


ij{x,Hi) = e{Ta,p,y,fJ.i)Ibb{Ta), i = (30) 

for the ROM. The rationale behind this choice is that in the arc simulations related to switchgear 
devices one can assume that initially the boundaries of the computational domain are at am¬ 
bient temperature and radiate the corresponding intensity into the volume filled with plasma. 
Additionally, the boundaries can not reach temperatures higher than about 3,000 K (depending 
on the material), such that the related radiated heat is anyway negligible with respect to the 
intensity emitted in the plasma volume where current is flowing. 


3.2.2 Model parametrization and model set 


The choice of the model parametrization, i.e. of functions /a and /g in (24), is crucially important 


for the accuracy of the obtained results and for the solution of the numerical optimization. For 
the first aspect, one shall choose a rich enough family of functions, such that the data from the 
FOM can be reproduced with small errors. For the second aspect, the best choice would be a 


parametrization leading to a convex model set (26), so that sequential quadratic programming 


and line search algorithms m can be used efficiently. For the equivalent absorption coefficients 
a, i.e. function fa, a choice that meets both requirements is a piecewise affine parametrization, 
where 6 a{p, V, Pi) is, for each band i = 1,..., M, a vector containing the values of the coefficients 
a at a finite number of pre-defined temperature nodes T^, k = such that T^-i < 

Tfc, k = 2 ,..., Nk, chosen by the user (e.g. equally spaced): 


da{p,y,Pi) = 


Oa{Ti,p,y,pi) 


^a{TNk,P,y,Pi) 


i = l,. 


,M. 


The temperature nodes must include the values at the boundaries 0 of the domain T, i.e. 

Tl = Tmin and Tm^ = Tmax- 


Then, for a given temperature value T £ T, the function fa in (|24|) is computed by interpolating 

(31) 


linearly among the values of 9a corresponding to the neighboring temperature nodes: 

fa{T,Oa{p,y,pi)) = \{T)'^eb{p,y,Pi), 


where 


A(r) = 


1 - 


T - Tk-i 

Tk — Tk-i 

T - Tk-i 

Tk — Tk-i 


E 


(32) 


0 

and Tk-i, Tk are two subsequent nodes such that T E [Tfc-i, Tk]. Using the parametrization 
(31), the number of optimization variables introduced in the problem is M Nk- The flexibility of 
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the approximating function can be increased by increasing the number of temperature nodes 
Nk- Moreover, note that fa in (31) is a convex combination [9] of the elements contained in the 
parameter vector. Therefore, with this parametrization, the constraints (23) can be enforced by 
imposing them just on the values at the nodes, i.e.: 


0 < a{Tk,p,y,^^i) < 1, A: = 1,... ,iVfc, i = 1,... ,M . (33) 

As regards the parametrization of the equivalent emissivity functions, e{T,p,y, pi), an addi¬ 
tional requirement is to retain a physical link between the emissivities of each band in the ROM 
and the original absorption spectrum, such that each pi accounts for the intensity contributed 
by a precise frequency interval defined by suitable frequency cuts. As already mentioned, pre¬ 
vious contributions in the literature actually consider such a partitioning of the EM spectrum 
Equh], however the choice of the frequency bands is not trivial and not systematic, so that 
one has to proceed with a trial-and-error approach. With the technique proposed here, one can 
optimize directly with respect to such frequency cuts, hence obtaining a systematic method for 
the band-averaging. This can be done by using the following model parametrization for the 
equivalent emissivities: 


fe(T, 9eip, y, Pi)) = jiT, Oeip, y, Pi)) = j e(r, v)dv, (34) 

where e(T, v) is the spectral emissivity ([^ and the optimization variables are the upper bound¬ 
aries (e.g. in Hz) of each frequency band: 


Oe{p,y,Pi) = 


9e^^{p,y,pi) 

m-i) 


G 


dM—1 


For the hrst and last bands, i.e. pi and pM, the equivalent emissivity is computed using (34) with 
/io —^ 0 and pM = +inf, respectively. In this way, the equivalent emissivity of each band is, by 
construction, equal to the one pertaining to the frequency interval y,/Xj_i), (p, y, pj)]. 


The constraints (22) defining the model set (i.e. the second and third rows in (26)) can be easily 
taken into account by a set of linear inequalities: 


Ghf{p,y,Pi) > 0 

Ghf{p,y,Pi+i) > 9hf{p,y,Pi), i = -2 ■ 


(35) 


The constraints (35) impose an increasing ordering of the boundaries. By construction, the 
resulting bands are not overlapping and they cover the whole spectrum: these features auto¬ 
matically enforce the constraints ( |22[ ). We note that with this parametrization, the number of 
free variables that define the band partition is equal to M — 1. Therefore, considering also the 
parameters pertaining to the equivalent absorption coefficients, the total number of optimization 
variables in the problem is equal to M + M — 1. About the total number of constraints, in 
virtue of (33) we have a set of 2M linear inequalities, while (35) amounts to M —1 additional 
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inequalities, for a total of 2 M Nj. + M — 1. These inequalities altogether define indeed a convex 
set (in particular a polytope) where the solution of the optimization problem is confined. 


About the computation of 7 (T, 0e(p, y, /i*)) in (34), note that this function can be also written 


as: 


where 


l{T,0e{p,y,yi)) = fj{T,ee^^{p,y,fii)) - f^{T,ee^^{p,y, pi-i)) 


(36) 


MT,e)= / e{T,v)du. 


(37) 


The function M'*' —>■ [0,1] can be conveniently pre-computed and stored, so that the 

computation of 'y{T,9^{p,y, pi)) can be very efficient. 


3.2.3 NLP solution and computational aspects 

After choosing the temperature profiles and the model parametrization as described in sections 


3.2.1 and 3.2.2, respectively, it can be shown that the cost function J (28) is twice differen¬ 


tiable. Thus, considering also the convexity of the model set, the optimization problem 
can be solved for a local minimum with constrained sequential quadratic programming (SQP) 
techniques m, like the one implemented in Matlab® function fmincon. The computational 
efficiency and stability of the numerical optimization depend on a number of aspects, which are 
briefly mentioned here. 

Computation of the cost function. SQP solvers employ (as most NLP solution algorithms) 
an iterative strategy where the cost function might need to be evaluated hundreds of times. It is 
thus imperative to speed-up the computation of such a function which, in our case, implies the 
computation of the L—dimensional vectors AIj{9) = Wj{Itot,j{Tj) — Itot,j{Tj, 6)), j = 1,..., Lt, 
see (27)-(28). The intensity profiles given by the FOM, Itot,j{Tj), can be compnted once and 
used in all the function evaluations, since they are not changing with 9. On the other hand, 
the intensity given by the ROM, Itot,j{Tj,9), has to be computed iteratively during the opti¬ 


mization. This implies evaluating the selected functions (24) and simulating the ROM model. 


As mentioned briefly in section 2.3, the computational times of such operations improve signif¬ 


icantly by using the discretized version of the ROM (18) and by considering the values of a 


instead of d, as optimization parameters, since one avoids the computation of the exponential in 


(20). The values of a can be then recovered from the optimal solution by inverting such equa¬ 


tion. Another approach to speed up the computations is to pre-compute and store the vectors 
X{Tj{xi)), i = 1,..., L, j = 1,..., Lt (32) corresponding to the employed temperature profiles, 


which are needed to compute the piecewise affine functions fa{T,9^{p,y,pi)) (31), since these 
vectors do not depend on the optimization variables but just on the temperature distributions 


Gradient computations. Another aspect that greatly influences the efficiency (and accuracy) 
of the optimization algorithm is the computation of the gradients of the cost and constraint 
functions. While the gradient of each constraint is trivial to compute, since only linear equalities 


and inequalities are present (as shown in section 3.2.2), the gradient of the cost is more difficult 
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to obtain, but we show here that it can still be derived analytically. First of all, consider that: 


Lx L 


j=l 1=1 


(38) 


hence 


VeJ(0) = 


dJ 


dJ 


Lx L 

El 

i=i i=i 


EE -2wj {hot, j{xi) Itot,jixi, 0)'j ^ 0 hot,j{xi, 0)■ (39) 


We thus focus on deriving a general expression for the gradient V 9 hot,jixi, 0) of the total radiated 
intensity at the /—th position node for the j—th temperature profile, which can then be employed 


considering the vector C in (13), we have: 


to compute the gradient of J using (39). Exploiting the second of the model equations (18) and 


M 


^ dhot,jixi, 6) — ^ ohjxij iJ-ii 0) 1 


(40) 


i=l 


where ij{xi, is the intensity pertaining to the band //j at the position step xi when the 


temperature profile Tj is considered. Moreover, by using the hrst of (18) together with (19) and 


(24) we obtain, for each I = 2,..., L: 


\/eIj{xuiXi,e) = VefaiTj{xi_i),eaip,y,iJ.i)) Ijixi_i,iJ,i,9) 

+ (1 - faiTjixi_i),9a{p,y, Pi))) y efe{Tj{xi-i),9^{p,y^pi)) hyiTj{xi_i)) 

- ^efa{Tj{xi-i),9a{p,y, Pi)) fe{Tj{xi_i),9e{p,y, Pi)) Ibb{Tjixi_i)) 

(41) 


The gradients Ve/a and Ve/g in (41) depend on the model parametrization. For the piecewise 
affine function /a, the computation is straightforward, e.g. (from (|31[)): 


'^efa{Tjixi-i),9a{p,y,Pi))) = Ve {X{Tj{xi-i))9a{p,y, pi)) = 


KTjixi-i)) 


0 


(42) 


where the zeros in vector Ve/a correspond to all the other parameters in vector 9 but the set 
9a(j),y,Pi) pertaining to the band. 

As regards function /g, whose parameters are the boundaries of the spectral bands, the com¬ 
putation of the related gradient Ve/g involves computing the derivative of function f^{T,9) in 


(37), which is readily obtained as df,y{T,9)/d9 = e{T,9) (^. More specihcally, considering 
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and (36)-(37) we have: 


'^efe{Tjixi_i),ee{p,y,i^i))) = V0-f{Tj{xi_i),ee{p,y,Pi)) 

0 


-e{Tj{xi-i), {p, y, Aij-i)) 
e{Tj{xi-i),ee^j {p,y, Pi)) 


0 


(43) 


where, in a way similar to (42), the zeros in vector Ve/g correspond to all the other parameters 
in vector 6 but the pair y, pi) pertaining to the band. For the first 

and last bands, i.e. i = 1 and i = M respectively, the gradient is computed by replacing the 
spectral emissivity pertaining to (p, y, yi_i) (resp. Oenj^iVi Pi)) with zero, since these two 
boundaries are not free variables, as commented above in Section [3.2.2 


The last ingredient needed for the gradient computation is the initialization of Vglj{xi, pi,9) 


conditions (30): 


for / = 1, which is required in (41); this is readily done by computing the derivative of the initial 


^ 0 lj{xi,Pi, 6 ) = V 0 fe{Ta, 6 e{p,y,Pi)) hb{Ta), i = 1,... ,M, j = 1,... ,Lt 


(44) 


By using (41)-(44), one can recursively compute the values of V 0 lj{xi, pi,9) for each position 
node xi, each temperature profile Tj and each band pi, and then use ([38|)-(|40[) to compute the 
gradient of the cost function. The recursive gradient computation can be done together with the 
simulation of the ROM, which is anyway needed to compute the cost function. The alternative 
to such an approach would be to estimate the gradient by means of numerical differentiation, e.g. 
finite difference approximation, which would imply computing + 1 times the cost function for 
a single gradient estimate. This is done by default in most available optimization routines if no 
gradient is provided, however it is subject to numerical errors and it requires significantly more 


time than the exact computation derived above. Therefore, the use of (38)-(44) to compute 


'V 0 J{ 9 ) is both more accurate and more efficient than numerical differentiation. 

Regularization and warm start. Other important issues that can arise in the numerical 


solution of (29) are related to the quality of the obtained optimizer and the stability of the 
optimization algorithm. One aspect that is relevant for the subsequent use of the derived 
absorption coefficients in CFD simulations is the smoothness of the obtained functions a and 
9 with respect to their input arguments, T,p,y,pi. About the temperature dependence, if 
piecewise affine functions are used one can penalize the variation of a (or equivalently d) directly 
in the optimization, by augmenting the cost function with a regularization term, as follows: 


J{9) = V AlJ(0)A/j(0) + /3A0^A0, 


i=i 
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Figure 9: Example of absorption coefficients and emissivities for a ROM with two bands for pure air. 


where 


M = 


Gd,2{P^y,fJ-i) - Od,i{p,y,yi) 
^d,ne^{P,y,yi) - Od,ng^-l{P,y,Pl) 
^d,2{p,y,PM) - 9d,i{p,y,PM) 

{p,y,PM) - {p, y, pm) 


and /3 is a weighting factor that can be used to tradeoff the variability of the absorption coef¬ 
ficients with respect to temperature with the goodness of the fitting criterion. The gradient of 
this modified cost can be computed with a straightforward extension of equation (39). Abont 
the smoothness with respect to p and y, recall that we assnmed that several ROMs are com¬ 
puted independently by gridding the V and y domains (see section 3.1). Thus, a possible way 
to influence the dependence of the coefficients on pressure and composition is to warm-start 
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Figure 10; Example of absorption coefficients and emissivities for a ROM with two bands for pure air. 
The gray lines correspond to the curves shown in Figure]^ 


the optimization, by initializing the optimization parameters with the solution obtained from 
neighboring values of p,y. This approach usually yields good results, as we show in the next 
section with some examples. 


4 Results 


We applied our method to compute the equivalent absorption coefficients and the partition of the 
EM spectrum in frequency bands for several gaseous media, ranging from pure air to mixtures of 
silver (or copper), air and hydrogen, and carbon dioxide and copper. As an example, Figures]^ 
show the results for pure air, using M = 2 bands in the ROM. For each temperature profile T, 
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the weight wj in (27) has been chosen by taking the average value of the corresponding intensity 
prohle computed with the FOM. From these results, it can be noted that for low pressure values 
(Figure]^, in a temperature range of about 12-1410^K the first of the two bands accounts for 
more than 80% of the total black-body intensity, however with an absorption coefficient of the 
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Figure 11: Optimal partitioning of the frequency spectrum for pure air and (a) two bands, (b) three 
bands. 


order of 10“^ m“^. Considering that the steady-state intensity is reached after about 3/a of 
propagation distance within the medium, this band can be considered to be transparent with 
respect to the space-scale of interest (i.e. about 10“^ m). On the other hand, in the same 
temperature range the second band has low emissivity (accounting for the remaining 20%), but 
large absorption coefficient, hence it is able to reach in short distance the corresponding fraction 
of black-body intensity. When large absorption coefficients are present (again relative to the 
considered space-scale), the band is said to be diffusive. In practice, however, the two bands 
emit (or absorb) often similar values of intensity in the same distance, since the transparent one 
contributes a small fraction of a very large intensity value, while the diffusive one contributes a 
large fraction of a small intensity value. 

The situation can change significantly for higher pressure values: in this case, the diffusive 
bands can have larger emissivity values at high temperature, meaning that they reach in short 
distance a significant portion of the black-body intensity (see Figure [To]). This implies, roughly 
speaking, that at these pressure and temperature values a larger quantity of radiated heat is 
redistributed to the surroundings and transferred to the walls. Similar considerations as the one 
just presented apply for all the other mixtures that have been considered. 


Figure 11-(a) shows the dependency of the optimal partition of the EM spectrum, i.e. of 


parameter (P) yj/^i) pressure. It can be noted that, as pressure increases, the boundary 
between the two bands shifts slightly and gradually towards lower frequencies, meaning that, for 
a given temperature, the second band accounts for larger fractions of the black-body intensity. 
This effect is clear also from the courses of the emissivities in Figures [9 10 In Fig. (HKb), 
the optimal partition achieved with M = 3 bands is shown. It can be noted that the variation 
of the frequency cuts with pressure is less marked in this case, in particular for the boundary 
between the first and second bands. An example of the corresponding absorption coefficients 
and emissivities is shown in Fig. 


12 


To give an idea of the accuracy achieved by the ROMs computed with the proposed method 
with respect to the FOM, we present two further examples. The first one is related to the 
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Figure 12: Example of absorption coefficients and emissivities for a ROM with three bands for pure air. 
Pressure values: 2.5110® Pa (solid lines), 3.98 10® Pa (dashed lines), 6.3110® Pa (dotted lines), 110® Pa 
(dash-dotted lines). 


radiated heat intensity along a single line of propagation, with the temperature profile shown 
in Figure 13 )a), through a mixture of 25% air, 50% copper and 25% hydrogen at 2.5 10^ Pa. 


Such a profile is a realistic temperature distribution that can take place in the plasma generated 
during the switching process of a circuit breaker, and it was not part of the data used to derive 
the ROM coefficients. Two ROMs are considered, one with two bands and the other with three 
bands. The resulting distributions of the radiated heat intensity are shown in Figure [T^b). R 
can be noted that both ROMs are able to reproduce quite accurately the intensity profile given 
by the FOM, with the third-order one being slightly more accurate. This example thus shows 
that increasing the number of bands gives in general higher accuracy, but the gain is smaller and 
smaller (in line with the considerations of section 2.2) and usually it is not worth using more 
than three bands, due to the increased computational load in the CFD simulations. 

The second example is related to the use of the derived ROM for a discrete-ordinate method 
(DOM, see e.g. [26]) simulation, where we want to compute the radiated heat intensity that 
reaches a wall in front of a column of hot gas composed by 75% copper and 25% air, at 10® Pa. 
Such a setup is described in Figure [l4)(a)-(b) . 

This example is more meaningful for the sake of CFD simulations of a plasma in real devices, 
where not just the intensity along a single line but the net total intensity obtained by integrating 
over all directions has to be computed, for each one of the finite volumes or elements that 
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Figure 13: (a) Temperature profile used to compare the radiated heat predicted by the FOM with that 
predicted by the ROM. (b) Comparison between the radiated heat intensities predicted by the FOM and 
the ones given by two ROMs with order 2 and order 3, respectively, using the temperature distribution 
of Figure (a). Mixture: 25% air, 50% copper and 25% hydrogen at 2.5 10® Pa. 

(a) (b) 




Figure 14: Setup for a simulation of the radiated heat intensity using the discrete ordinate method. A 
column of hot gas (75% copper and 25% air at 10® Pa) is standing in front of a wall and we want to 
compute the radiated heat that reaches the wall, (a) Scheme of the setup, the wall of interest is parallel 
to the (Y, Z) and contains the point (1.5 10“^, 0, 0). The column is represented by blue circles and the 
position of its center in {X,Y, Z) is shown as a solid black line, (b) Course of the temperature along a 
line parallel to the (X, Y) plane, passing through the center of the column, the highest temperature value 
corresponding to the column center. 


partition the computational domain. The intensity distribution (in W/m^) on the wall is shown 
in Figure 15 It can be noted that the ROM is able to capture well the qualitative behavior 
given by the FOM, with an average error of about 10%. Such a value is indeed very small as 
compared with the accuracy that can be achieved with other approaches, e.g. using constant 
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Figure 15; Simulation results for the hot column example with the DOM method. Radiated heat 
intensity (W/m^) that reaches the wall computed with the full-order model (left) and with a reduced- 
order model (right) with three bands and piecewise affine parametrization with 30 temperature nodes in 
the range [300, 410"^] K. 


absorption coefficients over frequency bands. 

Finally, as regards the computational load, here are some indications obtained by running 
the order reduction algorithm on a workstation equipped with 12 Intel Xeon® cores at 2.4 GHz 
each and 52GB of RAM in total. The average time required to compute one set of coefficients 
for a given pressure and composition pair {p,y) was 0.23 hours for order two models, and 0.32 
hours for order three models, in both cases with a temperature discretization with 50 nodes 
for the absorption coefficient, i.e. 101 and 152 optimization variables, respectively (see section 


5 Conclusions and future developments 

We presented a new approach to derive reduced order models of the radiated heat intensity 
through a participating gaseous medium. The approach is based on nonlinear identification 
and numerical optimization methods, and it can deliver models with low complexity but still 
highly accurate with respect to the original full absorption spectrum of the considered mixture. 
We also introduced a system-theoretic perspective of the phenomenon, as well as a frequency- 
domain analysis that justifies the use of low-order models to predict the total radiated intensity. 
We presented several comparisons between the full-order model and the reduced-order ones, 
obtained with our approach, and discussed the required computational effort. With the proposed 
technique, one can reliably and systematically derive radiation models that can be used in GFD 
simulations of a plasma, since both the average absorption coefficients and the frequency bands 
are a result of the numerical optimization. The accuracy of the model depends mainly on the 
accuracy of the base data that the model approximates. Future development efforts shall thus 
be devoted to assess and improve such base data. 
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